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Abstract 

A pressure-based multi-block computational method is developed for solving the 
incompressible Navier-Stokes equations in general curvilinear grid systems. The scheme is based 
on the semi-implicit type flow solver with the staggered grid. Issues concerning the mass and 
Tnn mMitnm flux treatments at the discontinuous grid interface are addressed. Systematic 
numerical experiments for different interface treatments involving (i) straightforward 
interpolation, (ii) globally conservative scheme, and (m) locally conservative scheme have been 
conducted. It is demonstrated that mass conservation has to be maintained locally, at the grid 
interface, with accuracy compatible with that of the scheme used in interior domain. Direct 
interpolation or globally conservative interface treatment of mass flux can not yield solutions 
with desirable accuracy. 



1. Introduction 

The numerical solution of fluid flow equations generally requires the generation of a grid 
for the region of interest. For many engineering problems with complex geometries, the 
generation of a single structured grid to cover the whole domain with the desirable grid 
distribution can be very difficult. For some flow problems with multiple length scales, it is also 
hard to generate a single structure grid to resolve all the flow features with reasonable grid 
points. These difficulties can be overcome to a limited extent by applying sophisticated grid 
generation schemes to construct a single grid with suitable characteristics. However, the degree 
of satisfaction achievable with such a process is highly problem dependent. To simplify this 
problem, it is becoming more common to use several grids at once, each with a reg ular grid 
structure. The various grids may either overlap in an irregular fashion or patch together. Farh 
component grid covers a relatively simpler geometry and can be generated individually. Since 
grid lines need not be continuous across grid interfaces, local grid refinement and adaptive 
redistribution can be conducted more easily to improve the solution accuracy. In the region with 
high flow gradients, by using such flexible grid layouts, the total storage and CPU time can be 
reduced while achieving the desirable solution accuracy. 

In order to make good use of the multi-block method, one need to appropriately handle 
the grid interface treatment associated with the flow solver. Because the grid lines may not be 
continuous across the block interfaces, interpolation and data communication methods have to 
be devised to transfer the information between blocks. These methods should be preferably easy 
to implement while maintaining good efficiency and desirable accuracy. Besides these 
requirements, for many fluid flow problems containing varying flow gradients, it is often 
important to use the conservative interface procedure to ensure that the physical laws are 
satisfied there[l,2]. Clearly, these considerations impose serious constraints on the construction 
of an interface scheme; furthermore, in some cases the above requirements can conflict with one 
another. Simultaneous achievement of both conservation and accuracy can be a very difficult 
task. 


Some progress has been made in this area, with different goals, for different fluid 
physics and multi-block arrangements, including patched and overlapped grids. Patched grids 
are individual grid blocks of which two neighboring blocks are joined together at a common grid 
line without overlap. With overlapped grids the grid blocks are partially superimposed on each 
other to cover the region of interest. Patched grid is relatively easy for data structure 
management, but it still has some difficulty with grid generation because of the interface 
constraints. Overlapped grid has more flexibility with grid generation, but, in general, its data 
management is more complex. For both grids arrangements, issues of interface treatment 
regarding both conservation laws and spatial accuracy need to be addressed. In the compressible 
flow regime, Rai [3,4] has developed conservative interface schemes for Euler equations 
calculation on the patched grid in a general curvilinear coordinate framework, for both explicit 
and implicit time integration schemes. Chesshire and Henshaw [5] have developed an overlapped 
grid generation method and a set of data structure, and solved the compressible Navier-Stokes 
equations. They treat the grid interfaces using interpolation without fluxes conservation. Meakin 
[6] has investigated the spatial and temporal accuracy of overlapped grid methods for invicid 
moving body problems, using tri-linear interpolation for grid interface treatment. He suggests 
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that grid resolution is the primary issue. Attempts have been made by Moon and Liou [7] and 
Wang and Yang [8] to devise conservative interface schemes for overlapped grids. However, 
the issues of the importance of conservative interface treatment versus solution accuracy, and 
the requirements in local and global conservation have yet been clearly addressed. 

In the incompressible flow regime, some progress has also been achieved. Hinatsu and 
Ferziger [9] have solved the unsteady Navier-Stokes equations in complex geometry using an 
explicit method. They does not enforce the mass conservation at the grid interface. Yung et al. 
[10] have solved the Navier-Stokes equations with a semi-implicit algorithm in the patched 
curvilinear grid system. However, the fluxes in general are not treated conservatively across grid 
interface. Lai et al. [11] have solved Navier-Stokes equations with the patched curvilinear grid 
system. They treats all the diffusion terms implicitly across the grid interfaces by defining 
individual nodal points for each block in the common areas. This way, while easy to implement, 
does not guarantee the satisfaction of the conservation laws across interface either. Meakin and 
Street [12] have solved the unsteady Navier-Stokes equations using the overlapped grid and 
applied the method to three-dimensional environmental flow problem. Tu and Fuchs [13] have 
developed a computational methodology for Navier-Stokes with emphasis on using an overlapped 
grid technique and multigrid method, and applied the method to the unsteady three-dimensional 
internal engine flow simulation. Neither work enforces conservative interface treatment. 
Henshaw [14] has used a fourth-order accurate method to solve incompressible Navier-Stokes 
equations on overlapped grids. In his method, the discrete divergence of the velocity field is not 
exactly zero, and a damping term is needed to stabilize the computations. Wright and Shyy [15] 
developed a pressure-based multi-block method for solving the incompressible Navier-Stokes 
equations on domains composed of an arbitrary number of overlapped grid blocks in the 
Cartesian grid system. A locally conservative internal boundary scheme, with first order 
accuracy, is devised to ensure that global conservation of mass and momentum fluxes are 
maintain ed. This methodology has also been extended to the curvilinear coordinates [16]. 
Clearly, although much progress has been achieved in the composite grid with complex fluid 
flow problems computed, some important and fundamental issues still need to be investigated 
further. For example, for incompressible flow problems, the mass conservation should be strictly 
maintaine d over the entire flow boundary; otherwise, the compatibility condition will not be 
satisfied. For the multi-block method, with which the discontinuous grid interface can be 
introduced into the flow domain, the conservative treatment of the mass flux along the grid 
interface may be critical for obtaining a converged solution. However, a interface scheme of 
lower order accuracy than the interior nodes, although conservative, may not be desirable from 
the resolution viewpoint, either. These issues have not been thoroughly investigated. 
Furthermore, for incompressible flow problems, the pressure may be only known up to an 
arbitrary constant. When the solver is applied to the multi-block grid, the pressure in different 
grid blocks may be independent. The coupling of the pressure at the interface can affect the 
solution process. 

In the present study, the issues of conservation and accuracy of the interface treatment 
in a multi-block method is investigated for incompressible viscous flow computations on general 
curvilinear staggered grids. Effects of non-conservative and conservative treatments, either only 
globally along the whole interface or locally for each computational cell, for mass flux across 
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the grid interface on solution accuracy are focused on. The methods are tested for a recirculating 
problem with multi-block grids. A variety of numerical experiments with different grid 
resolutions are conducted, and the results are compared with the benchmark solution to assess 
the issues involved. The momentum flux treatment involving convection, pressure and viscous 
terms at the interface for momentum equations are also discussed. 

2. Governing equations and numerical algorithm 

The governing equations adopted in this study are the 2-D steady state, incompressible, 
constant property, Navier-Stokes equations along with the corresponding form of the continuity 
equation. 


4<p«> + -r<pv> - o 

ox dy 


(la) 


-i(p ««) * -i(puv) - -&■ ♦ <■ 

dx dy dx dx dx dy dy 


(lb) 


— (puv) + — (pw) = + — (p— ) + —([!—) 

dx 8y By ax dx dy dy 
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With the introduction of coordinate transformation £ = £(x,y), ij = y(x,y), the equations above 
are cast in the curvilinear coordinates, 


^-(plD ♦ x-(pP) - 0 

as an 


(2a) 


^(pVu) + ±( (2b) 

(2c) 


where 

U = uy, -vx,, V = vxj - uy t 

qi = x 2 , +f„ q> = x t x, + y*y„ % = x\ + y 2 t 

J = x t y, - x,y t 

A staggered grid system is adopted, as shown in Fig. 1. The scalar variables are located 
at the center of the four adjacent grids. Both u and U are located at the midpoint of the east and 
west faces of the control volume. Both v and V are located at the midpoint of the north and 
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south faces of the control volume. In terms of the notation shown in Fig. 1, for a node p 
enclosed in its cell and surrounded by its neighbors n, s, e, w, the finite-difference 
approximation to the conservation laws can be performed by taking the integral of momentum 
equations over the control volume and discretizing it. By arbitrarily taking A£=Atj = 1, the 
resulting u- and v-momentum equations yield 

[p^ + Vj^“r¥,)Hw + [pFh-y 5 P-j(-^ 2 u 5 + 9 3 “ n )] I" - 0 (3a) 


[pUv-x^-^jiq^-q^)] |t + [pKv+x^--t(-g 2 v { +tf 3 v n )] |* - 0 ( 3b ) 

The equations above can be put into a general difference form: 

\‘ w + [p^-yC-^^^t,)] I* = S J (4) 

where <t> is the general dependent variable and S is the source term. With appropriate finite 
Hi fferftnr.p. schemes representing the convective and diffusive terms at the control volume 
boundaries, the discretized equation relating the variable at a central point p and its neighboring 
values is obtained [17], 


a p$ p = (5) 

where a’s are the coefficients resulting from the numerical schemes chosen in course of 
discretization. The pressure term and cross-derivative portion of the viscous terms due to the 
non-orthogonal coordinates are treated here as the source term S*. The continuity and momentum 

equations can be used to formulate a pressure correction equation. The pressure correction, p' , 
is used to update the pressure field, and in conjunction with a velocity correction formula to 
obtain a velocity field satisfying the continuity equation at convergence. The discretized form 
of the pressure correction equation is presented as: 

a pP'p " a ePi +a „Pl+ a J>'s a sPs +S r (6) 




a +a +a +a 

t w * s 


Sp = (p tO w ~(p u*) e + (p ^),-(p p*), 

where the superscript * designates the intermediate solution, and p' designates the correction 
made. A detailed discussion of this algorithm can be found in Refs. [17,18] and will not be 
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repeated here. In single-block grid computations, the solution procedure is as follows: (a) The 
momentum equations are solved to obtain the Cartesian velocity components with the given 
pressure field. When solving the momentum equations, the contravariant velocity components 
U and V are calculated after updating each of the velocity components, (b) With die updated U 

and V, the pressure correction equation is solved to obtain p' . (c) The velocity and pressure 
fields are updated with the solution of the pressure correction equation. Procedures (a)-(c) are 
repeated until the momentum and continuity equations are simultaneously satisfied to the required 
degree of accuracy. The procedure above is a semi-implicit procedure, which is reduced to the 
SIMPLE procedure if the Cartesian grid is adopted [19]. With curvilinear coordinates, a 
combined use of both Cartesian and contravariant velocity components is devised. In particular, 
while the Cartesian components are the primary dependent variables in momentum equations, 
the contravariant components are corrected first in the pressure correction step to ensure the 
satisfaction of mass continuity. Because of the semi-implicit nature of the procedure, when the 
algorithm is applied to the multi-block grid computation, different computational strategies 
between different grid blocks can be adopted: (1) Within each grid block, procedures (a)-(c) are 
repeated several times to update the solutions without updating the boundary and interface. 
Afterwards, the computations are conducted in the neighboring blocks. Such a block to block 
cycling procedure continues until convergence is achieved in all blocks. (2) While solving each 
of the u-, v- momentum and pressure correction equations, iteration is conducted from block to 
block, without updating other equations. In other words, procedures (a)-(c) are conducted in the 
outer loop with die multi-block computation embedded within each differential equation. In this 
study, strategy (1) is adopted, since it allows the interface treatments among different equations 
to be handled together. The computation continues until the mass and momentum residuals, 
normalized by the mass and momentum fluxes entering the whole domain, in each block meet 
the convergence criterion. More information regarding both the composite grid techniques and 
the pressure-based algorithm discussed here cjm be found in Shyy [20]. 

3. Grid interface treatment 

The basic arrangement in the present multi-block grid system is that there is at least one 
grid layer overlap, for each block, between the adjacent blocks. However, The grid lines from 
both blocks need not be continuous in the overlapped region. At grid interfaces, the boundary 
conditions for each block are needed in order to solve the equations. Because the grid interfaces 
are not the physical boundaries, the interface boundary values are not known a priori and must 
be obtained as a part of the whole solution. Some interpolations are needed to acquire the 
intermediate interface boundary conditions between neighboring blocks. Then, a block to block 
iteration is conducted to obtain the solution in the whole domain. 

3.1 Interface treatment for the momentum equations 

For the momentum equations, e.g., the u-momentum equation, the fluxes 

= ptf“ + (7a) 
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F\ v = fiVu-yp- 


(7b) 


could be interpolated directly from neighboring blocks at the interfaces. Since they include 
velocity gradient and metric terms, such a direct interpolation does not results in a conservative 
solution. In the present treatment, direct interpolation of dependent variables, with adjustment 
of pressure to maintain the global momentum conservation at grid interfaces, is adopted. The 
quadratic or the linear interpolation can be used [5]. With the staggered grid arrangement, the 
interface treatments for the east and west boundaries are different from those for the north and 
south boundaries. Here, the east boundary is defined as the grid boundary with grid index 
i = W’ the west boundary as with index i= 1 . The north and south boundaries are the other 
sides with j =j max and j = 1 respectively. 

First consider the evaluation of the flux E at the east or the west interface. Due to the 
utilization of the staggered grid, the control volume face at which the u-momentum flux is 
evaluated is not at the grid boundary, as illustrated in Fig. 2. The geometric quantities q„ and 
J are well defined within the current block. The velocity component u at the grid interface 
boundary is interpolated from the neighboring block to evaluate the viscous term in E. The 
pressure term at the control volume face is also available within the current block. 

Second, consider the evaluation of the flux F at the north or south interface. The control 
volume face coincides with the grid interface. The geometric quantities and the viscous terms 
are not well defined within the block. With present treatment, the interface grid lines are 
extended from the original block to intersect with the neighboring grid lines to form fictitious 
interpolation points as illustrated in Fig. 3. Then the geometric quantities are obtained based on 
linear interpolation. The u components are interpolated at the fictitious points and both 
convection and viscous terms can be evaluated from the values obtained from the original block 
and the fictitious nodes. Generally, for curvilinear grid, the pressure term will appear in the flux 
F if y { is not zero. Because the pressure is not defined at the control volume face, either it can 
be extrapolated from inside the current block or interpolated from the adjacent block. The latter 
one is adopted in this study. With the present solver, the boundary condition for the pressure 
correction equation is of Neumman type, which means the pressure is known up to an arbitrary 
constant. In the present approach, a pressure adjusting method is used to connect the pressure 
filed of different blocks consistently, which is based on the total momentum flux balance across 
the interface [15]. This adjustment is conducted between the block to block iterations to keep 
the pressure fields compatible between the blocks. The interface treatment above also can be 
applied to the v-momentum equation similarly. It is noted that with the present pressure 
adjustment procedure, the total momentum fluxes across the grid interface are uniquely 
determined for both blocks, resulting in a conservative treatment. 

3.2 Interface treatment for the mass flux 

In general, two types of boundary conditions can be used for the pressure correction 
equation. If the pressure is known at the boundary, the value of the pressure correction there is, 
of course, zero. If, instead of the pressure itself, the velocity component normal to the boundary 
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is prescribed (such as the inflow and solid wall conditions), then, due to the staggered grid 
arrangement, the Neumman type of boundary condition is imposed on the pressure correction. 
The value of pressure correction at the boundary is no longer needed. Consider the pressure 
correction Eq. (6), suppose the normal velocity component is known at the south boundary, the 
equation for a control volume next to that boundary becomes 

a P Pp = o t p' t *a w p , w *a ll p , m +S p 
a p = a a +a„+a n 

S p = (ptO w -(pt/V(pK),-(pn. 

Where (pF), is the mass flux from the south boundary and is known, U* and V* are the 

intermediate values to be further corrected by the pressure correction p 1 . As stated above, p t 
does not appear because no "correction" is necessary on the south boundary. Since the continuity 
equation is not solved explicitly, the outflow boundary condition should always satisfy the global 
mass conservation. Usually, a global mass correction is applied to outflow boundary during the 
solving process to help maintaining the global mass conservation. This correction will not affect 
the final solution when the process is finally converged. Otherwise, the process will either 
converge very slowly or even fail to converge [21]. Now consider the case that the south 
boundary is the grid interface. In this situation, (pP) x is not known as a priori during the 
solution process. It can be obtained from the intermediate solution in the overlapped region of 
the neighboring block at the previous iteration. But with the discontinuous grid interface, direct 
interpolations generally will not satisfy the mass conservation across the grid interface. Non- 
conservative error will be introduced at the grid interface and the magnitude of error will depend 
on the order of the interpolation method used at the grid interface. In this study, the linear and 
quadratic non-conservative interpolations are implemented to investigated the effect of non- 
conservative error from the grid interface on the solution. Also, several conservative treatments 
of different local and global natures, and of different orders of accuracy are implemented and 
are tested. 

(a) Interpolation with global correction: 

First, a non-conservative interpolation is used to obtained the mass flux at the interface 
from the neighboring block. The total flux evaluated from neighboring block is denoted as F„ 
the total interpolated flux is denoted as F c . Second, the mass deficit A F = F a -F e is computed. 
Then the mass flux at each current interface control volume face is added with AF/N (suppose 
there are N control volumes on the interface) so that the total flux across the interface is 
conserved globally. But this conservation is not enforced locally, and its effect on the solution 
in the whole domain need to be assessed. 


(8a) 

(8b) 

(8c) 
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(b) Piece-wise constant interpolation: Such a scheme is used in Refs. [3,15], which by nature 
is locally conservative, but with only first-order accuracy. 

(c) Linear (or quadratic) interpolation of mass flux with local conservative correction: 

Suppose a portion of an interface, corresponds to the width of a single control volume 
of the coarse grid 1, and to the width of several control volumes of the fine grid 2, indexed from 

i=l to imax. This scenario is shown schematically in Fig. 4, where V e is the contravariant 
velocity component normal to the grid face, normalized by the control volume face length S c , 

at the coarse grid control volume face, and represents the contravariant velocity component, 
no rmalize d by the control volume face length 5^, at the fine grid control volume face. From 

fine grid to coarse grid, V e can be obtained as 


^ imax __ 

V, - E V/. (9) 

i-l 

Equation (9) is for the constant density flow. It can be easily extended to account for density 
variation. For the sake of simplicity, this aspect is neglected in the current discussion. In this 
way, the flux into the coarse grid is uniquely determined from the corresponding fine grid fluxes 

conservatively. Conversely, given the coarse grid flux V c S c , the conservation constraint yields 
the follows: 


imax 


E % - v . s . 


( 10 ) 


In this situation, conservation does not provide unique values for the V fi . A certain distribution 

has to be chosen to decide each . As a first approximation, is obtained using the linear 
(or quadratic) interpolation of the normalized contravariant velocity in the coarse grid. The 

interpolated value is denoted as . With the first approximation, Eq.(9) is not satisfied. Then 

the fine grid fluxes are scaled so that the total flux obtained is V C S C . Accordingly, the values 
at the fine grid boundary are computed as follows: 


v s fi = - V.S. 

fi fi imax 

E \^\ s n 


1*1 


( 11 ) 


From Eq.(ll), it can be seen that Eq.(9) expressing flux conservation from coarse grid to fine 
grid is satisfied and the flux distribution is close to that determined by linear (or quadratic) 
interpolation. 
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If the grid interface is not exactly matched in the coarse to fine grid manner, as 
illustrated in Fig. 5, then, a split-merge procedure is devised in the above mentioned spirit, as 
shown in Fig. 6. The control volume in grid 2 can be split into smaller subcontrol volumes. In 
this manner, from control volumes in grid 1 to subcontrol volumes in grid 2, the locally 
conservative treatment can still be applied on a cell-by-cell basis, after the interpolation has 
beat conducted. After that, the fluxes at the split control volume faces will be merged back to 
get the flux at the original control volume face. From grid 2 to grid 1, this same split-merge 
treatment is also applicable. The only extra work in this procedure is to create some arrays to 
store the intermediate information. Thus, the interface can be treated no matter what kind of 
interface arrangement is encountered. This linear (quadratic) interpolation with local correction 
treatment is not limited to the mass flux, it also can be applied to the momentum flux 
conservative treatment. 

3.3 Data structure 

Because of the possible grid discontinuity, a set of data structure has to be devised to deal 
with information transfer between different grid blocks. The data structure should not only 
provide the necessary information for the chosen interface treatment but also be as simple as 
possible. In the present method, all the data structure information is based on the grid 
distribution. The minimum requirement for the grid interface between the two adjacent grid 
blocks is that there exists at least one grid layer overlap. For each block, each side of the grid 
boundaries may be divided into several segments corresponding to different boundary condition 
types. For each segment, the indices of die starting and the ending grid points, the boundary 
condition type and the corresponding values of the dependent variables on the physical 
boundaries are assigned. If the segment is a grid interface, the block number of the adjacent 
block, the identification of overlap direction (in £ or q direction), and the number of overlap 
layers are assigned as well. Based on the input above, a series of intersection tests are conducted 
to provide the information for the coefficients for interpolation and conservative treatment. These 
coefficients are stored for later use. 


4. Numerical experiments and discussions 

The test problem is the lid-driven cavity flow problem with Re= 1000. First, a 3-block 
discontinuous Cartesian grid configuration with different grid resolutions is used. A grid system 
of 41x21, 81x13 and 41x11 grid points for block 1, 2 and 3, respectively, as shown in Fig. 7 
is first employed. Here, the interface between block 1 and 2 coincides with the cavity horizontal 
center line. The second grid system doubles the grid resolutions of the first system, and has 
81x41, 161x23 and 161x21 grid points in each block, respectively. The third grid system doubles 
the grid resolution in x direction of the second grid system and has 161x41, 321x23 and 161x21 
grid points in three blocks. To clarify the terminology, the composite grid for the whole flow 
domain is called a "grid system" here. Each grid system consists of several blocks. The first grid 
system described above is denoted the coarse grid system, the second the median grid system, 
and the third the fine grid system. The three grid systems share the same topological 
characteristics in each block. We have created these three systems to investigate the interplay 
of interface treatment and overall grid resolutions. The grid layouts are not ideally suitable for 
the present recirculating flow, they are purposely set up to test the relative merits of different 
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interface treatments. For all the test cases presented in this work, the second-order central 
difference is used for convection, diffusion and pressure terms. 

Case 1. The Cartesian velocity components u and v and the contravariant velocity 
components U and V are linearly interpolated at the grid interfaces. Total momentum fluxes 
across the interface are conserved via the pressure adjustment. No mass conservation is enforced 
across the grid interfaces. The computations are conducted over the coarse, median and fine grid 
systems. For all three grid systems, the normalized residuals for u and v momenta are below 
10 4 . The normalized mass residuals reach down to 1.5X10* 3 , 3.3xl0 4 and 1.3xl0 4 for the 
coarse, median and fine grid systems, respectively, and stabilize at those levels. Figures 8a and 
8b show the u -component distributions along the vertical center line and v-component 
distributions at the horizontal center line; they are compared with the corresponding benchmark 
solution reported by Ghia et al.[22]. It can be seen that the solutions for all the three grid 
systems have substantial discrepancies with respect to the benchmark values. Although in 
general, the solutions improve as the grids are refined, the overall performance of all three 
systems are unsatisfactory. It is noted that for this problem computed with single grid distributed 
uniformly, a 81x81 grid system can yield very accurate solution already [23]. Accordingly, it 
is unsatisfactory to observe that both median and fine grid systems, even with resolutions better 
than the 81x81 single uniform grid, still do not yield accurate solutions. Figure 8c shows the 
interface pressure distributions computed on block 1 and 2. The pressure in each block has been 
adjusted according to the total momentum flux balance at the interface. The absolute value of 
the pressure has no practical meaning. For a well converged solution, the pressures from 
different blocks are expected to have the same distributions at the interface; in the current case, 
even with the fine grid system, the discrepancy between the two pressure distributions is 
obvious. 

Case 2. Both the Cartesian velocity components u and v and the contravariant velocities 
U and V are quadratically interpolated, but still without enforcing mass conservation. This 
procedure is of higher interpolation accuracy than case 1; however, both are without 
conservative treatment. The motivation here is to test the role played by the inteipolation 
accuracy, and to investigate that for viscous flow, whether conservative treatments are necessary 
for obtaining satisfactory solutions. For all three grid systems, the residuals for the momentum 
equations are below 10 4 . The mass residuals stabilize at the level of 2.3xl0‘ 3 , 4.3X10 4 and 
7.5xl0" 5 for the coarse, median and fine grid systems, respectively. Figures 9a and 9b show the 
u-component and v-component distributions respectively for the three grid systems. Overall, 
these solutions are more accurate than those shown in Fig. 8, and the solution on the fine grid 
system shows obvious improvement over the coarse and median grid systems. However, 
discrepancies still exist between the present and the benchmark solutions. Figure 9c shows the 
pressure distributions from block 1 and 2 at the interface. The discrepancy between the two 
distributions also still exist. From cases 1 and 2, it can be seen that the solution can be improved 
with the grid refinement. The overall solution accuracy is better with quadratic interface 
interpolation than with linear interpolation, which is expected because the quadratic interpolation 
increases the order of interpolation accuracy and should reduce the non-conservative error for 
mass flux across the grid interfaces. But even with the fine grid system, the solution of quadratic 
interpolation is still not satisfactory. Since with both interpolation schemes, the solutions are not 
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as accurate as the single continuous grid results, the source of this inaccuracy must come from 
the non-conservative interface treatments across discontinuous grid blocks. 

Case 3. Since the non-conservative linear and quadratic interpolations for the 
contravariant velocities cannot lead to a satisfactory solution even for a very fine grid, a 
conservative interface treatment is then tested. In this case, the Cartesian velocity components 
u and v are linearly interpolated. The contravariant velocities U and V are linearly interpolated 
first, followed by a global correction procedure for mass flux as discussed previously. The 
computation is conducted for the coarse grid system only. The momentum and mass flux 
residuals reach to the level of Id 5 . The results are shown in Figs. 10a, 10b and 10c. Obviously, 
tiie solution is not satisfactory. It appears that the conservative interface treatment conducted at 
the global level does not improve the solution accuracy. 

Case 4. The Cartesian velocity components u and v are still linearly interpolated. The 
contravariant velocities U and V are interpolated based on the piecewise constant formula, which 
by nature is locally conservative with first-order accuracy. This treatment is implemented to 
investigate the effect of the local conservation on the solution. The computations are conducted 
for the three grid systems. For the coarse and median grid systems, the residuals for the 
momentum and mass fluxes reach down to Id 5 . But for the fine grid, the solution does not 
converge. Figures 11a and lib present the u and v component distributions for the coarse and 
median grids. Both u and v profiles agree well with the benchmark solutions. Figures 11c and 
lid exhibit the pressure distributions at the interfaces of block 1 and 2 for two grid systems. The 
pressures in the interface region obtained on different blocks follow each other generally well, 
except that high wave number oscillations appear on block 2. The cause of this nonphysical 
oscillation comes from the fact that the mass flux distributions in the overlapping region of the 
fine grid cells are assigned according to the piecewise constant formula. This distribution 
formula results in a series of stair-step mass flux profile on the fine grid block, forcing the 
pressure field to oscillate in response to the non-smooth mass flux distribution. This same reason 
is also probably responsible for the nonconvergence of the fine grid system. 

Case 5. The Cartesian velocity components u and v are linearly interpolated. The 
contravariant velocities U and V are linearly interpolated first, followed by a local correction 
to maintain the cell-by-cell mass conservation across the interfaces. The residuals for momentum 
and mass fluxes reach down to 10 -5 . Figures 12a and 12b show the u and v component 
distributions for three grid systems. The solution for the coarse grid system shows a very small 
discrepancy compared to the benchmark solution; the solutions on the median and fine grid 
systems agree very well with the benchmark solution. Figure 12c displays the pressure 
distributions at interface of block 1 and 2 for the fine grid system. The pressure distributions 
from the two adjacent grid blocks conform to each other very well. Clearly, local conservation 
holds a key to produce satisfactory solution accuracy. 

Case 6. The Cartesian velocity components u and v are quadratically interpolated. The 
contravariant velocities U and V are quadratically interpolated first, and then a local correction 
is applied to maintain the mass conservation across the interfaces. Again, the computations are 
conducted over the three grid system and the residuals are at the level of Id 5 . Figures 13a and 
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13b show the u and v-component distributions. Good agreements, comparable to Case 5 (Fig. 
12) between the current solutions and the benchmark solution are observed. Figure 13c presents 
the pressure distributions at block 1 and 2 interface for the fine grid system. The pressure 
distributions at the interface also conform to each other very well. For the present flow problem 
both the linear and quadratic interpolation with local correction give the satisfactory solutions. 
Since the discretization scheme for the interior cells are second-order accurate, it appears that 
a linear interpolation (aided with the follow-up local conservation treatment) procedure is 
sufficient. 

Case 7. Finally, the same flow problem is investigated with a 3-block curvilinear grid, 
which has 81x41, 161x23 and 81x21 grid points for block 1 ,2 and 3, respectively, as shown 
in Fig. 14. The interface treatment is the same as that in Case 5, and the momentum and mass 
residuals reach to lO" 5 . Figure 15 demonstrates the u velocity component distribution at the 
vertical center line and compares it to the corresponding benchmark solution. Good agreement 
is obtained; illustrating the flexibility of the curvilinear grid system. 

5. Conclusions 

A pressure-based multi-block computational method has been developed for solving the 
incompressible Navier-Stokes equations in a general curvilinear grid system. For the momentum 
equations, the pressure fields between two adjacent blocks allow an arbitrary jump, which can 
be adjusted by conserving the total momentum fluxes across the block interface. The importance 
of main taining local mass flux conservation across the interface with certain accuracy is 
illustrated through a series of numerical experiments. Specifically, both the linear and quadratic 
interpolations without conservative measure for the mass flux at the interface cannot lead to the 
desired solution even for very fine grid. Linear interpolation with global correction for mass flux 
does not improve the solution, either. The piecewise constant treatment can improve the solution 
accuracy due to its conservative nature, but creates artificial pressure oscillation due to non- 
smooth mass flux distribution in the fine grid block. Nevertheless, it illustrates the importance 
of local mass flux conservation across the interface. Both linear and quadratic interpolations, 
with local conservative correction prove to be good choices. The fact that their solutions are of 
very comparable accuracy indicates that (1) the interface treatment needs not to be of higher 
formal order of accuracy than the interior schemes, and (2) the local conservative mass flux 
treatment at interface holds a key for composite grid computations. 
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(i) Configuration of a staggered grid system. 



(ii) Finite difference grid representstion : 

(a) physical plane; (b) transformed plane. 


Fig. 1 Staggered grid and notation for both 
physical and computational domains 
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block 1 


south inteface- 
for block 1 
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u 
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block 2 


Fig. 3 The ficticious points arrangement for u at the south interface of block 1 . 
The dashed lines are extended from the solid lines in block 1 and intersect with 
a grid line in block 2. These intersection points will be used to estimate the flux 
for block 1 . 



Fig. 4 Notations used in two-block interface (subscript c and f designate coarse 
and fine grid quantities, respectively) 
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grid 1 


grid 2 


Fig. 5 General interface configuration: no distinct coarse or fine grid. 



grid 1 


grid 2 


Fig. 6 Split control volume configuration: control volumes in grid 2 are split 
into subcontrol volumes 
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Fig. 7 Three-block grid with 41x21, 81x13 and 41x11 nodes 
for block 1 ,2 and 3, respectively. 
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Fig. 8c The pressure distribution at the interface between 
block 1 and 2 on the fine grid system. 

Fig. 8 Solution profiles based on linear interpolation for U and V 
without mass conservation correction. 
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Fig. 9a The u velocity along the vertical center line. 



Fig. 9b The v velocity along the horizontal center line. 
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Fig. 9c The pressure distribution at the interface between 
block 1 and 2 on the fine grid system. 


Fig. 9 Solution profiles based on quadratic interpolation for U an V 
without mass conservation correction. 
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Fig. 10a The u vwlocity along the vertical center line 



Fig. 10b The v velocity along the horizontal center line. 
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Fig. 10c The pressure distribution at the interface between 
block 1 and 2 on the coarse grid system. 


Fig. 10 Solution profiles based on linear interpolation with global correction 
of mass conservation for U and V on the coarse grid system. 
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Fig. 11a The u velocity along the vertical center line. 



Fig. 1 lb The v velocity along the horizontal center line. 
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Fig. 1 lc The pressure distribution at the interface between 
block 1 and 2 on the coarse grid system. 



Fig. lid The pressure distribution at the interface between 
block 1 and 2 on the median grid system. 

Fig. 11 Solution profiles based on piecewise constant interpolation for U and V. 
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Fig. 12a The u velocity along the vertical center line. 



Fig. 12b The v velocity along the horizontal center line. 
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Fig. 12c The pressure distribution at the interface between 
block 1 and 2 on the fine grid system. 


Fig. 12 Solution profiles based on linear interpolation with local correction 
of mass conservation for U and V. 
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Fig. 13c The pressure distribution at the interface between 
block 1 and 2 on the fine grid system. 


Fig. 13 Solution profiles based on quadratic interpolation with local correction 
of mas conservation for U and V 
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block 1 



Fig. 14 Three-block curvilinear grid with 81x41, 161x23 and 81x21 
nodes for block 1, 2 and 3, respectively. 
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Fig. 15 The u velocity along the vertical center line. Linear interpolation 
with local correction of mass conservation for U and V. 
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